Long gravitational-wave transients and associated detection strategies for a network of 

terrestrial interferometers 
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Searches for gravitational waves (GWs) traditionally focus on persistent sources (e.g., pulsars or 
the stochastic background) or on transients sources (e.g., compact binary inspirals or core-collapse 
supernovae), which last for timescales of milliseconds to seconds. We explore the possibility of long 
GW transients with unknown waveforms lasting from many seconds to weeks. We propose a novel 
analysis technique to bridge the gap between short 0(s) "burst" analyses and persistent stochastic 
analyses. Our technique utilizes frequency-time maps of GW strain cross-power between two spa- 
tially separated terrestrial GW detectors. The application of our cross-power statistic to searches 
for GW transients is framed as a pattern recognition problem, and we discuss several pattern- 
recognition techniques. We demonstrate these techniques by recovering simulated GW signals in 
simulated detector noise. We also recover environmental noise artifacts, thereby demonstrating a 
novel technique for the identification of such artifacts in GW interferometers. We compare the 
efficiency of this framework to other techniques such as matched filtering. 

PACS numbers: 95.55.Ym 



I. INTRODUCTION 

Historically, searches for gravitational- wave (GW) 
transients fall into one of two categories: searches for 
"bursts" whose precise waveforms we cannot predict and 
searches for compact binary coalescences, whose wave- 
forms can be predicted (at least for the inspiral part). 
Typically, burst searches focus on events with < 1 s du- 
rations and, indeed, there are many compelling models 
for short GW transients (see [l[ and references therein). 

In this paper, we put the spotlight on long GW tran- 
sients whose durations may range from many seconds to 
weeks. Astrophysical GW emission scenarios for long 
transients exist (e.g., 2 5]), but their characteristics 
have not previously been broadly addressed and no data- 
analysis strategy has been proposed for such events until 
now. (In addition to this work, see recent developments 
in Q.) Most of the GW emission models we consider are 
burst-like in the sense that the signal evolution cannot be 
precisely predicted, however, we refer to them as "tran- 
sients" to avoid connoting that they are short-duration. 

In Sec. [ill we survey a range of mechanisms for GW 
emission that may lead to long transients. These in- 
clude long-lived turbulent convection in protoncutron 
stars (PNSs), rotational instabilities in rapidly spinning 
PNSs and in double neutron-star merger remnants, mag- 
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netoturbulence and gravitational instabilities in gamma- 
ray burst (GRB) accretion torii, r-modes associated with 
accreting and newborn neutron stars, as well as, per- 
haps more speculatively, pulsar glitches and soft-gamma- 
repeater (SGR) outbursts. 

In Sec. IIII1 we introduce an analysis framework utiliz- 
ing frequency-time (/i)-maps of GW strain cross-power 
created using data from two or more spatially separated 
detectors. The framework is extended to include multi- 
ple detectors, and we show that it is a generalization of 
the GW radiometer algorithm Q- In Sec. HV] we com- 
pare /i-cross-power maps of GW data (time-shifted to 
remove astrophysical content) with Monte Carlo simula- 
tions of idealized detector noise. We shall see that GW 
interferometer data is well-behaved enough that thresh- 
olds for candidate events can be estimated analytically 
(in at least one case). 

In Sec. [V] we use /i-cross-power maps to cast the 
search for long GW transients as a pattern recognition 
problem. For the sake of concreteness, we consider two 
algorithms: a "box search" Q and a Radon algorithm |{|. 
In Sec. lVIl we demonstrate the Radon alg orithm (as well 
as the "locust" and Hough algorithms [lfj]) to identify en- 
vironmental noise artifacts in LIGO environmental mon- 
itoring channels — a novel technique for the identification 
of such artifacts in GW interferometers. In Sec. I VIII we 
describe how our framework is related to other detection 
strategies such as matched filtering. Concluding remarks 
are given in Sec. I Villi 
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model 



PNS convection core collapse 

rot. instabilities BNS coalescence, core collapse, isolated NS 

r-modes core collapse, isolated NS 

disk instabilities BNS coalescence, core collapse 

high-e BH binaries BBH coalescence 

pulsar glitches isolated NS 

SGR flares isolated NS 

TABLE I: Models of long GW transients with associated 
sources. BNS and BBH stand for "binary neutron star" and 
"binary black hole" respectively. 



II. ASTROPHYSICAL SOURCES OF LONG GW 
TRANSIENTS 

In this section, we review a variety of emission mech- 
anisms for long GW transients. Most of the mechanisms 
we consider (summarized in Tab. [XJ) are associated with 
one or more of three types of objects: core-collapse su- 
pernovae (CCSNe), compact binary inspirals, or isolated 
neutron stars. 



A. Core-collapse supernovae and long gamma-ray 
bursts 

There is tremendous electromagnetic observational ev- 
idence connecting both CCSNe and long gamma-ray 
bursts (GRBs) to the core-collapse death of massive stars 
(see, e.g., [U). Both are ultimately powered by the re- 
lease of gravitational energy, but the precise mechanism 
by which gravitational energy is converted into energy 
of ejecta and radiation is uncertain in both phenomena 
(see, e.g., (lT| - [l3j and references therein). However, all 
modern models of CCSN and long-GRB central engines 
involve violent non-spherical dynamics, making both sys- 
tems prodigious emitters of GWs. 

The GW signature of CCSNe (recently reviewed in Q) 
may be composed of contributions from rotating col- 
lapse and core bounce II 4J] , post-bounce protoneutron 
star (PNS) convection (2l.ll5l.fl6j. neutrino-driven convec- 
tion and the standing-accretion-shock instability (SASI) 
[l7l - [l9j , PNS pulsations [2(| , nonaxisymmetric rotational 
instabilities (both dynamical and secular) [2l|, [22[ , asym- 
metric neutrino emission [2| 
[23|-[22} , magnetic stresses 



in rotating PNSs (see, e.g., 




aspherical outflows 
and r-modc pulsations 
Depending on the 
particular CCSN mechanism in operation, some emission 
processes may dominate while others are suppressed [l3j . 

Currently, there are two favored scenarios for the long- 
GRB central engine. In the collapsar scenario [30j . mas- 
sive stars collapse to black holes either without an initial 
CCSN exploding or via fallback accretion after a success- 
ful, but weak explosion. The millisecond-protomagnctar 
scenario (3ll.l32j relies on highly magnetized, nascent neu- 



tron stars. In both cases, any long-GRB activity is pre- 
ceded by stellar collapse and a post-bounce phase during 
which a PNS exists and GW emission occurs in very simi- 
lar fashion to regular CCSNe. In the collapsar scenario, a 
black hole with an accretion disk forms. Magnetohydro- 
dynamical processes and/or neutrino pair annihilation 
powered by accretion and/or by the extraction of black 
hole spin energy eventually launch the GRB jet. GWs 
may be emitted by disk turbulence and disk instabilities 
that may lead to clumping or disk fragmentation (HQ . In 
the millisecond-magnetar scenario, a successful magneto- 
rotational CCSN explosion (see, e.g., [H, [H|) occurs, 
after which a high-Lorentz-factor outflow is driven by 
the millisecond protomagnetar. GWs may be emitted by 
convectivc/mcridional currents and dynamical and secu- 
lar nonaxisymmetric rotational instabilities in the proto- 
magnetar @, Hl|. 

In CCSNe and in the CCSN-phasc of long GRBs most 
GW emission processes last until the onset of the CCSN 
explosion or until PNS-collapse to a black hole, and hence 
they have a short duration of order < 1 — 2 s [HI, H?} ■ 
Exceptions arc PNS convection, secular rotational insta- 
bilities including r-modes and long-GRB disk/torus in- 
stabilities. We discuss these below and provide order-of- 
magnitudc estimates of their emission characteristics in 
the time and frequency domains. 



1. Protoneutron star convection 

If a CCSN explosion occurs, a stable PNS is left be- 
hind and will cool on a Kelvin-Helmholtz timescale (see, 
e.g., [II])- Fallback accretion [39|, or, perhaps, a late- 
time hadron-quark phase transition (e.g., (!o|) may lead 
to PNS collapse and black hole formation. If the PNS 
survives, a powerful convective engine, driven by ther- 
mal and lepton gradients may continue to operate for 
possibly tens of seconds in the cooling phase [H, |41H43| , 
making it a long GW transient source. 

GW emission estimates from PNS convection are based 
on results of simulations that track only the early phase 
(< 1 s after core bounce) 0, HH, [44[ , yet they have found 
a number of robust features that translate to later times. 
PNS convection occurs at moderate to high Reynolds 
numbers, hence, is turbulent and leads to an incoher- 
ent, virtually stochastic GW signal. Its polarization is 
random in the nonrotating or slowly rotating case, but 
may assume specific polarization due to axisymmetric 
rotationally-driven meridional currents in rapidly spin- 
ning PNSs (an effect that remains to be studied in com- 
putational models). In the phase covered by current 
models, typical GW strains are h ~ 3x 10~ 23 at a galactic 
distance of lOkpc @, EB]- ("Strain" refers to the strain 
measured at Earth; strain amplitude scales like the in- 
verse of the distance from the source.) 

While on short timcscales, the GW signal of PNS 
convection will appear almost as a white-noise burst, 
its time-frequency structure is non-trivial, exhibiting a 
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broad spectral peak at 0(100 Hz), which shifts to higher 
frequencies over the course of the first second after core 
bounce @, [Hj] . This chirp- like trend is likely to continue 
for seconds afterward as the PNS becomes more com- 
pact. It should be nearly independent of the size of the 
convectively unstable region, since the eddy size will be 
set by the local pressure scale height (see, e.g., (45j). 

Based on the ~ 1.2 s evolution of a PNS model by 
[l|| Hl|, we expect a total emitted energy of ~ 1.6 x 
10 °MqCT. Assuming that convective GW emission 
continues with comparable vigor for tens of seconds, we 
can scale this to E GW ~4x 10~ 9 (At/30 s) M Q c 2 . 



2. Rotational instabilities 

Most massive stars (> 98% or so, see [46[ and refer- 
ences therein) are likely to be slow rotators, making PNSs 
with birth spin periods of ~ 10 — 100 ms. GRB progen- 
itors, however, are most likely rapidly spinning, leading 
to PNSs with birth spins of Oil ms) and rotational ki- 
netic energy of up to 10 52 erg [33j, enough to power a long 
GRB through protomagnetar spin-down as suggested by 
the protomagnetar model [3l|, [32[ . 

PNSs, like any self-gravitating (rotating or nonrotat- 
ing) fluid body, tend to evolve toward a state of minimal 
total energy. PNSs are most likely born with an inner 
core in solid-body rotation and an outer region that is 
strongly differentially rotating [46[. Magnetorotational 
instabilities (see, e.g., 1471 ') and/or hydrodynamic shear 
instabilities (see, e.g., [48] ) will act to redistribute an- 
gular momentum toward uniform rotation (the lowest- 
energy state). The latter type of instability may lead 
to significant, though short-term t < Is, nonaxisym- 
metric deformation of parts of the PNS and, as a con- 
sequence, to significant GW emission @, HHH^]. PNSs 
in near solid-body rotation that exceed certain values 
of the ratio of rotational kinetic to gravitational en- 
ergy, T/|W|, may deform from an axisymmetric shape 
to assume more energetically-favorable triaxial shape of 
lowest-order I = m = 2, corresponding to a spinning 
bar. Such a bar is a copious emitter of GWs. At 
T/|W] > 0.27, nonaxisymmetric deformation occurs dy- 
namically, but will not last longer than a few dynamical 
times of O(ms) (see, e.g., [5(1 HU) due to rapid angular 
momentum redistribution, and hence we do not consider 
its GW emission in this study. 

At T/| W| > 0.14, a secular gravitational-radiation re- 
action or viscosity-driven instability may set in, also lead- 
ing to nonaxisymmetric deformation. The timcscalc for 
this depends on the detailed PNS dynamics as well as 
the details and strength of the viscosity in the PNS. It is 
estimated to be 0(1 s) for both driving agents, but the 
expectation is that gravitational-radiation reaction dom- 
inates over viscosity [H, HH . The secular instability has 
the potential of lasting for ~ 10 — 100 s [H, HH , and hence 
it is of particular interest for our present study. 

Once the gravitational-radiation reaction instability 



sets in, the initially axisymmetric PNS slowly deforms 
into I = m = 2 bar shape and, in the ideal Dedekind 
ellipsoid limit, evolves toward zero pattern speed (angu- 
lar velocity = 0) with its remaining rotational energy 
being stored as motion of the fluid in highly noncircu- 
lar orbits inside the bar [22|, [52[ . GW emission occurs 
throughout the secular evolution with strain amplitudes 
h proportional to Q 2 and to the ellipticity e, charac- 
terizing the magnitude of the bar deformation, leading 
to an initial rise of the characteristic strain followed by 
slow decay as decreases [U HH ■ We expect charac- 
teristic strain amplitudes, defined as h c = tqw/gw' 1 , 
of h c = O(10~ 23 - 10~ 22 Hz" 1 / 2 ) for a source located at 
100 Mpc and the emission is expected to last at that level 
for O(100 s) [HH. The emitted GWs will be elliptically 
polarized. 



3. R-modes 

i?-modes are quasi-toroidal oscillations that have the 
Coriolis force as their restoring force. It was shown 
in [28l . Ilij that r-modes in neutron stars are unsta- 
ble to growth at all rotation rates by gravitational- 
radiation reaction via the secular Chandrasekhar- 
Friedman-Schutz instability [55l IBH ] . R- modes emit (at 
lowest order) current-quadrupole GWs with /gw = 
4/3 (Qns/27t) and typical strain amplitudes h ~ 4.4 x 
10- 24 a(O N s/v^G^) 3 (20Mpc/ J D) |H3|, where fl NS is the 
NS angular velocity, D is the distance to the source 
and p is the mean neutron star density. The parameter 
a € [0, 1] is the dimensionless saturation amplitude of the 
r-modes and its true value has been the topic of much 
debate. Most recent work suggests (see 0, [58{ and ref- 
erences therein) that a <C 0.1 and, perhaps, does not ex- 
ceed ~ 10~ 5 due to non- linear mode coupling effects [59| . 
Generally, r-modes are expected to be a source of very- 
long-lasting quasi-continuous GW emission, though long 
GW transients may be possible in the case of high satu- 
ration amplitudes (e.g., [60|V 

Potential astrophysical sources of GWs from r-modes 
are accre ting neutron stars in low-mass X-ray binaries 
(e.g., [58l. l6ll. l62l|) and, more relevant in the present con- 
text, newborn, rapidly spinning neutron stars [28l . l29l . l57l ]. 
In the latter, r-modes may play an important role in the 
early spin evolution [29|, |63|] . 

4- Accretion disk instabilities 

In the long GRB collapsar scenario, the central en- 
gine consists of a black hole surrounded by an accre- 
tion disk/torus [ill, H9|- The inner part of the disk 
is likely to be sufficiently hot to be neutrino cooled 
and thin [64| while the outer regions with radius r > 
50i?s = 100 GMbrc~ 2 are cooled inefficiently and form 
a thick accretion torus [H, [64|. A variety of (magneto)- 
hydrodynamic instabilities may occur in the disk/torus 
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leading to the emission of GWs. 

Piro and Pfahl [3] considered gravitational instability 
of the outer torus leading to fragmentation facilitated 
by efficient cooling through helium photodisintcgration. 
Multiple fragments may collapse to a single big dense 
fragment of up to ~ 1M© that then travels inward either 
by means of effective viscosity and/or GW emission. In 
both cases, the inspiral will last O(10— 100 s), making it 
a prime candidate source for a long GW transient. Piro 
and Pfahl predict maximum dimensionless strain ampli- 
tudes \h\ ~ 2 x 10- 23 (/ G w/1000Hz) 2 / 3 for emission of 
a system with a fragment mass of IMq, a black hole 
mass of 8M Q , and a source distance of 100 Mpc. The 
frequency will slowly increase over the emission interval, 
making the emission quasi-periodic and, thus, increasing 
its detectability by increasing its characteristic strain h c 
up to O(xl0 -22 ) at /gw ~ 100 Hz and a source distance 
of 100 Mpc. 

Van Putten, in a scries of articles (see, e.g., [f| 
loll loll ] and references therein), has proposed an ex- 
treme "suspended-accretion" scenario in which the cen- 
tral black hole and the accretion torus are dynamically 
linked by strong magnetic fields. In this picture, black 
hole spin-down drives both the actual GRB central en- 
gine and strong magnetoturbulence in the torus, leading 
to a time- varying mass quadrupole moment and, thus, to 
the emission of GWs. Van Putten postulated, based on a 
simple energy argument, branching ratios of emitted GW 
energy to electromagnetic energy of £gw/-Bem ^ 100 
and thus, -Eqw ~ few x 10 53 erg. These numbers are 
perhaps unlikely to obtain in nature, but the overall con- 
cept of driven magnetoturbulence is worth considering. 

In van Putten's theory, the magnetoturbulcnt torus ex- 
citations produce narrowband elliptically polarized GWs 
with a frequency between (1 ~ 2kHz)(l + z) for a GRB 
at a redshift of z [651 ] . The frequency is predicted to 
vary with time such that df/dt = const (6f|. GW emis- 
sion continues for approximately the same duration as 
the electromagnetic emission, lasting typically for sec- 
onds to minutes [ll[. With van Putten's energetics, a 
long GRB at a distance of 100 Mpc is predicted to pro- 
duce a strain of h ~ 10~ 23 , which is comparable to the 
expected Advanced LIGO noise at 1000 Hz J65|. Inte- 
grating many seconds of data, such a loud signal should 
stand out above the Advanced LIGO noise, making it 
likely that a strong statement can be made about this 
model in the advanced detector era. 



B. 



Postmerger evolution of double neutron-star 
coalescence 



In Subsec. Ill Al we discussed a variety of scenarios for 
long GW transients in the context of PNS and black- 
hole - accretion-disk systems left in the wake of CCSNe 
and in collapsars. A similar situation is likely to arise 
in the postmerger stage of double neutron-star coales- 
cence. The initial remnant will be a hot supermassive 



neutron star that, depending on the mass of the binary 
constituents and on the stiffness of the nuclear equation 
of state, may survive for hundreds of milliseconds (e.g., 
[67| and references therein). In these systems, many of 
the GW emission mechanisms discussed in the stellar col- 
lapse scenario may be active. Hence, it may be fruitful 
to search for long GW transients following observed in- 
spiral events as well as following short GRBs, (which are 
expected to be associated with binary inspirals). 

Inspiral events, however, need not invoke the forma- 
tion of a PNS in order to produce a long GW transient. 
Highly eccentric black hole binary inspirals are expected 
to produce complicated waveforms that are difficult to 
model with matched filtering and may persist for hun- 
dreds of seconds [H, [69| , and thus they are suitable can- 
didates for long GW transient searches. According to 
some models [70[ , a significant fraction of BBH form dy- 
namically with high eccentricities (e > 0.9) leading to 
an Advanced LIGO event rate of ~ 1 — 100yr _1 . Given 
such high rates, it is highly likely that these models can 
be thoroughly probed in the advanced detector era. 



C. Isolated neutron stars 

Isolated neutron stars are another potential source of 
long GW transients. In the following, we discuss pul- 
sar glitches and soft-gamma repeater flares as potential 
sources of GWs. 

Pulsar glitches are sudden speed-ups in the rotation of 
pulsing neutron stars observed by radio and X-ray obser- 
vatories. The fractional change in rotational frequency 
ranges from 10 -10 < A/// < 5 x 10 -6 , corresponding 



to rotational energy changes of < 10 erg [7l|, |72J. The 
speed-up, which takes place in < 2 min, is followed by a 
period of relaxation (typically weeks) during which the 
pulsar slows to its pre-glitch frequency [73| . 

The mechanism by which pulsar glitches occur is a 
matter of ongoing research [74| - l78j . and the extent to 
which they emit GWs is unknown. We therefore follow 
Andersson et al. [7^] and assume that the emitted en- 
ergy in GWs is comparable to the change in rotational 
energy. Given these energetics, and assuming a sim- 
ple exponentially-decaying damped waveform, a nearby 
[d = 1 kpc) glitch can produce, e.g., a O(10 s) quadrupole 
excitation with a strain of h ~ 8 x 10~ 24 at 3.8kHz [79i ]. 
This is about six times below the Advanced LIGO noise 
floor, which effectively rules out the possibility of detec- 
tion. A long GW event measured by Advanced LIGO 
and coincident with a pulsar glitch would therefore sug- 
gest a radically different glitch mechanism than the one 
considered in |79j . 

Flares from soft-gamma repeaters (SGRs) and anoma- 
lous X-ray pulsars ( AXPs) , which may be caused by seis- 
mic events in the crusts of magnetars, have also been 
proposed as sources of GWs. Recent searches by LIGO 
have set limits on lowest-order quadrupole ringdowns in 
SGR storms (U and in singlc-SGR events |l|. SGR 



5 



giant flares are associated with huge amounts of electro- 
magnetic energy (10 44 — 10 46 erg), and they are followed 
by a 0(100 s)-long tail characterized by quasi-periodic 
oscillations (see, e.g., [HJ). It is hypothesized that quasi- 
periodic oscillations following SGR flares may emit GWs 
through the excitation of torsional modes [83(. 

Current models of GW from SGRs/AXPs HMH are 
very preliminary, but even if we assume that only 0.1% of 
the 10 46 erg of electromagnetic energy in a nearby SGR 
flare is converted into GWs, then SGRs are nonetheless 
attractive targets in the advanced detector era. Current 
experiments have set limits on the emission of GW energy 
from SGRs (in the form of short bursts) at a level of 
Egw < 10 45 erg 0, Hl| (« 10% of the electromagnetic 
energy of a giant flare depending on the waveform type) . 
Since sensitivity to -Egw h 2 , it is likely that we can 
probe interesting energy scales (£"gw ~ 0.1% £em) in 
the advanced detector era. 



III. AN EXCESS CROSS-POWER STATISTIC 
A. Definitions and conventions 

Our present goal is to develop a statistic yr which can 
be used to estimate the GW power Hr (or the related 
quantities of GW fluence Ft and energy Er) associated 
with a long GW transient event confined to some set of 
discrete frequencies and times T. In order to define GW 
power, we first note the general form of a point-source 
GW field in the transverse-traceless gauge: 



h ab (t,x) = ^2 dfe£ b (n)~h A (f)e 



2irif(t+n-x/c) 



(3.1) 



Here O is the direction to the source, A is the polar- 
ization state and {e^ b } are the GW polarization tensors 
with Cartesian indices ab, (see App. IA II for additional 
details). Since Eg. 13.11 describes an astrophysical source, 
the Fourier transform of the strain h(f) is defined in the 
continuum limit. 

We now, however, consider a discrete measurement on 
the interval between t and t + T measured with a sam- 
pling frequency of f s , which corresponds to N s = f s T 
independent measurements. We utilize a discrete Fourier 
transformations denoted with tildes: 



» =W E fee" 3 "**/*- 

S 71 = 

N.-l 
k=0 



N.-l 



(3.2) 



The variable t — separated from other arguments with a 
semicolon — refers to the segment start time, as opposed 
to individual sampling times, denoted by t with no semi- 
colon, e.g., s(i). Variables associated with discrete mea- 
surements are summarized in Tab. [TTJ 



f s sampling frequency 

St = 1/ fs sampling time 

t segment start time 

T segment duration 

Sf frequency resolution 

N s number of sampling points in one segment 
TABLE II: Variables describing discrete measurements. 



The GW strain power spectrum (measured between t 
and t + T with a sampling frequency f s in a frequency 
band between / and / + Sf) is: 



H AA ,(t;f) = 2(h* A (t;f)h A ,(t;f)). 



(3.3) 



The factor of 2 comes from the fact that H AA > (t; /) is the 
one-sided power spectrum. 

It is convenient to characterize the source with a single 
spectrum that includes contributions from both + and x 
polarizations. We therefore define 



H(t;f) = Tv[H AA >(t;f)} 



(3.4) 



so as to be invariant under change of polarization bases. 
This definition is a generalization of the one-sided power 
spectrum for unpolarized sources found in 0, HH, |93[ . 

Our estimator ^Vr(^) utilizes frequency-time (ft)- 
maps: arrays of pixels each with a duration determined 
by the length of a data segment T and by the frequency 
resolution Sf. In Subsec. IIIID1 we describe how J^r can 
be constructed by combining clusters of /t-map pixels. 
We thereby extend the stochastic-search formalism de- 
veloped in 0, [92], |93| beyond models of persistent un- 
polarized sources to include polarized and unpolarized 
transient sources. In doing so, we endeavor to bridge the 
gap between searches for short 0(s) signals and stochas- 
tic searches for persistent GWs. We begin by considering 
just one pixel in the /t-map. 



B. A single /t-map pixel 



In App. IA 21 we derive the form of an estimator Y 
for GW power H(t; f) in a single /i-pixel by cross- 
correlating the strain time series si(t) and sj(t) from 
two spatially separa ted detectors, / and J, for a source 
at a sky position Cl [l!8| . We find that 



Y(t;f,Q) = Re 



(3.5) 



Qij(t;f,n)Cu{t;f) . 

Here, Cu(t; f) is the one-sided cross-power spectrum 

C IJ (t;f) = 2S}(t;f)Sj(t;f). (3.6) 

Meanwhile, Qu(t; f,Cl) is a filter function, which de- 
pends, among other things, on the source direction and 
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polarization. For unpolarizcd sources (see App. 

1 



Qu{t-,f,n) 



where eu(t;Cl) £ [0, 1], the "pair efficiency," is 
eiJ (t;ft) = iX;F/(f;ft)i^(t;A). 



(3.7) 



(3.8) 



Here F^(t; ft) is the "antenna factor" for detector / and 
Axu = xi — xj is the difference in position vectors of 
detectors / and J; (sec App. IA Pair efficiency is de- 
fined such that a GW with power H will induce I J strain 
cross-power given by eijH. It is unity only in the case 
where both interferometers are optimally oriented so that 
the change in arm length is equal to the strain ampli- 
tude. For additional details (including a derivation of 
the pair efficiency for polarized sources) see Apps. IA 11 
!A~2landlA~5l 

The variance of Y is calculated in App. IA 31 Then 
in App. IA 41 we show that the following expression for 
a Y (t; f, ft) (motivated by analogy with stochastic analy- 
ses [92j ) is an estimator for the variance of Y, 



1.x 



& Y (t-j,n) = ^QuitJM^^frfiP^it-J), (3.9) 



jadj/ 



where Pj dl is the average one-sided auto-power spectrum 
in neighboring pixels, 



Pf\t;f) = 2\s I (t;ff 



(3.10) 



The overline denotes an average over neighboring pix- 
els fll9j . 

From Eqs. 13.91 and 13.51 we define the signal to noise 
ratio SNR(i; /, ft) for a single /i-map pixel: 



SNR(t; /, ft) = Y(t; /, n)/a Y (t; /, ft) 

)u(t;f,£l) dj(t;f) 



=Rc 



(3.11) 



It depends on the phase of Qu(t;f, ft), but not on the 
magnitude. Thus, a single /t-pixel taken by itself con- 
tains no information about the polarization properties 
of the source, since the polarization does not affect the 
phase of Q. This degeneracy is broken when we combine 
/t-pixels from different times or from different detector 
pairs. 



C. Energy, fluence and power 

One of the most interesting intrinsic properties of a 
transient source of GWs is the total energy emitted in 
gravitational radiation, £?gw- By measuring -Eqw (and, 



when possible, comparing it to the observed electromag- 
netic energy, -Eem), we can make and test hypotheses 
about the total energy associated with the event as well 
as constrain models of GW production. Thus, it is use- 
ful to relate Y(t; f, ft) to Eqw and the related quantity 
of fluence. If the GW energy is emitted isotropically (in 
general it is not) then |94| . 



E GW = AnR' 



16ttG 



dt[hUt) + hi(t) 



(3.12) 



where R is the distance to the source. It follows that the 
equivalent isotropic energy is related to our cross-power 
estimator as follows: 



E GW {t;f,Q) = 47ri? 2 



AG 



{Tf 2 )Y(t;f,Q). (3.13) 



Ecw(t', /, ft) may contain significant uncertainty about 
the distance to the source or the isotropy of the GW 
emission. It is therefore useful to define a statistic that 
contains only uncertainty associated with the strain mea- 
surement. The natural solution is to construct a statistic 
for GW fluence, FQ W (t; /, ft), which is given by 



F GW (i;/,ft) 



~-Tf 



AttR 2 

7TC 



3 \ 



(3.14) 



4G )F (t ;/,ft). 



In the subsequent section, we show how multiple pix- 
els can be combined to calculate the average power in- 
side some set of pixels. The same calculation can be 
straightforwardly extended to calculate the total fluence. 
This is done by reweighting Y(t; f, ft) and a(t; f, ft) by 
(ttc 3 /4G)(T/ 2 ). Also Eqs. EHU and must be scaled 
by the number of pixels in a set, N; (otherwise we obtain 
average fluence instead of total fluence) . 



D. Multi-pixel statistic 

We now generalize from our single-pixel statistic to 
accommodate transients persisting over N pixels in some 
set of pixels, F. We define Hr to be the average power 
inside T, 



(3.15) 



A minimum-variance estimator for the GW power in T 
can be straightforwardly constructed from a weighted 
sum of Y(t; /, ft) for each pixel in V, 



3Mft) 



E ti/er y(t;/,n)ffy(t;/,n)- 



(3.16) 



Here we assume that the power is either evenly or ran- 
domly distributed inside T, which is to say (H(t; /)) = 



7 



(H(t; /')> = H and so (H r ) = H . Thus, 



The associated uncertainty is 



E t . /Gr Y(t;/,ft)<7y(t;/,ft)- 



-1/2 



\ E t;/e p^(i;/^)- 2 

E t;/e r(^( f ;/."))^( f ;/.")' 
E t;/e r^(^;/^)- 2 



(3.17) 



An 



<#r>. 



Here we have additionally assumed that there are no 
correlations between Y(t; /, Q) in different pixels. If 
the GW signal in different pixels is correlated, then the 
{Y(t;f, Cl)} are correlated and Eq. 13.161 should, in the- 
ory, be modified to include covariances between different 
pixels. In practice, however, the covariance matrix is not 
known, and so we must settle for this approximation, 
which gives the estimator a higher variance than could 
be achieved if the covariance matrix was known. 

The associated estimator for the uncertainty is 



-1/2 



(3.18) 



The choice of the set of pixels T to include in the sum in 
Eq. l3.16l is determined by the signal model. For example, 
a slowly varying narrowband signal can be modeled as a 
line of pixels on the ft-m&p. We explore this and other 
choices for T in greater detail in Sec. |V] 

The SNR for given a set of pixels V is given by 



SNR r (ft) 



y r (ft) 



(3.19) 



Since SNRp is the weighted sum of many independent 
measurements, we expect, due to the central limit theo- 
rem, that the distribution of SNRr will be increasingly 
well-approximated by a normal distribution as the vol- 
ume of r increases and more pixels are included in the 
sum fl~20l. 



E. Multi-detector statistic 

It is straightforward to generalize 5V for a detector 
network Af consisting of n > 2 spatially separated detec- 
tors. First, we generate n{n— l)/2 /f-maps for each pair 
of interferometers. Then we extend the sum over pixels 
in Eq. 13.161 to include a sum over unique detector pairs 
p(I,J): 



By construction, the expectation value is 



Hr 



(3.20) 



(3.21) 



<^(ft) = ( E E 
<p(i,j) ft 



(3.22) 



Adding new detectors to the network improves the 
statistic by mitigating degeneracies in sky direction and 
polarization parameters and also by improving sensitiv- 
ity to Hr by increasing the number of pixels contributing 

to 3^. 

F. Relationship to the GW radiometer 

The multi-pixel statistic is straightforwardly re- 
lated to the GW radiometer technique, which has been 
used to look for GWs from neutron stars in low-mass 
X-ray binaries Q. By constructing a rectangular set of 
pixels consisting of one or more frequency bins and last- 
ing the entire duration of a science run, we recover the 
radiometer statistic as a special case. 

It is instructive to compare the unpolarizcd radiometer 
statistic Q with our 



y rad (i; /,n) 



IIJ 



(*;/,«) 



x 



~s*At;f)Sj(t;f) 
7 (f;/,ft)il(/) 



(3.23) 



(3.24) 



Pi(f)Pj(f) 
e u (t;Cl)e 2mfCl - ASlj/c . (3.25) 



Here 7/ ,/(£;/, ft) is the so-called overlap reduction fac- 
tor, A t is a normalization factor and eij(t; £1) is the pair 
efficiency, which we define in Eq. |A"20l and Eq. 

There are two things worth noting here. First, the ex- 
tra factor of H(f)/ P](f)Pj(f) in the expression for Qjj 
does not appear in our expression for Qjj (see Eq. IA19I) . 
The factor of l/Pj(f)Pj(f) is proportional to o(/)~ 2 , 
and so it is analogous to the weighting factors in Eq. 13.161 
The difference is that f rad builds this weighting into the 
filter function whereas we opt to carry out the weighting 
when combining pixels. The factor of H(f) in Q/j d is 
the expected source power spectrum. When we choose 
a set of pixels T, we effectively define H(f) such that 
H(f) = const inside T and H(f) = outside T. 

Second, we note that apparently Qff oc e whereas 
our filter scales like Qjj oc 1/e. It turns out that both 
filters scale like 1/e because the radiometer normaliza- 
tion factor A oc e~ 2 . The historical reason for this is 
that the radiometer analysis was developed by analogy 
with isotropic analyses [921 ] . which includes an integral 
over all sky directions. The inclusion of 7(t; /, 0) in 
the expression for Qff serves to weight different direc- 
tions as more or less important just like the factor of 
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l/Pi(f)Pj(f) weights different frequencies. We opt to 
remove 7 in favor of e, which deemphasizes the analogy 
with the isotropic analysis in order to highlight the net- 
work sensitivity, which is characterized by e. 



G. Relation to other search frameworks 

This is not the first time that /i-maps of data have 
been proposed to search for GWs. The literature 
on this subject is extensive and diverse. We concen- 
trate on comparison with "excess power" methods, (see 
e.g., H, US llll ) ■ The key difference between our frame- 
work and others is that we cross-correlate data from two 
interferometers before they are rendered as_/i-maps. Pre- 
vious implementations such as (9F 



and 



instead form 



/t-maps by auto-corrclating data from each interferom- 
eter individually and then correlating regions of signifi- 
cance in these maps. For Gaussian noise, neither of these 
ways of combining data from different detectors is opti- 
mal. Instead, the optimal multi-detector method incor- 
porates both autocorrelated and cross-correlated com- 
ponents 0]. Real interferometric GW data, however, 
is not Gaussian. Rather, there is an underlying Gaus- 
sian component with frequent non-Gaussian bursts called 
"glitches." For situations of this type, our approach has 
two advantages. 

First, noise bursts in both detectors that coincide 
in time and frequency increase the false-alarm rate for 
statistics with auto-correlated components, but are sup- 
pressed in our cross-correlation analysis unless the wave- 
forms of the burst themselves arc correlated in phase like 
a true GW. Second, even when noise bursts are present, 
the pixel values in an /i-map of cross-correlated data are 
well approximated by a simple model. This is unlike ft- 
maps with auto-correlated components, for which there is 
a no simple description. Thus, while our statistic is sub- 
optimal for Gaussian data, we expect it to perform well 
for real interferometer data. Moreover, even in the case 
of Gaussian noise, we do not sacrifice much sensitivity 
compared to the optimal excess-power statistic, or even 
to matched filtering, as demonstrated in Sec. IVIII We 
compare the sensitivity of the cross-correlation statistic 
to other methods in Sec. IVIII 



IV. DISTRIBUTION OF SIGNAL AND 
BACKGROUND 

In order to determine if a candidate event warrants fur- 
ther examination, it is necessary to determine the thresh- 
old above which an event is elevated to a GW candidate. 
This threshold is usually phrased in terms of a false-alarm 
rate (FAR). In Sec. IIIII we argued that /i-maps of cross- 
power provide a convenient starting point for searches 
for long transients because cross-correlation yields a rea- 
sonably well-behaved SNR(i; /, Cl) statistic whose proba- 
bility density function (PDF) we can model numerically, 



thus allowing straightforward calculation of a nominal de- 
tection threshold in the presence of Gaussian noise. We 
now assess this claim quantitatively. 

We consider 52 s x 0.25 Hz pixels (created through a 
coarse-graining procedure described below), which are 
the intermediate data product in stochastic analyses 
such as LIGO's recent isotropic result [92} • There are, 
of course, other choices of pixel resolution, and differ- 
ent sources call for different resolutions. Typically one 
must balance concerns about the signal duration, the 
signal bandwidth and the stationarity of the detector 
noise. The PDF of SNR(t; /, Cl) for a single pixel cru- 
cially depends on details of the pixel size. E.g., the 
PDF of SNR(i; /, Cl) for coarse-grained pixels (described 
in App. IB II) is more nearly Gaussian-distributed since 
coarse-grained pixels are created by averaging over more 
than one frequency. Our goal here, therefore, is not an ex- 
haustive treatment. Rather, we aim to assess the agree- 
ment of data with our model using one pixel size, and in 
doing so, demonstrate how this assessment can be carried 
out in general. 

Our noise model assumes Gaussian strain noise, uncor- 
rected between detectors I and J, 



2(n I (t;f)nj(t;f))=6 IJ N I (t;f), 



(4.1) 



where Nj(t;f) is the one-sided noise power spectrum 
and hi{t; f) is the discrete Fourier transform of the noise 
strain time series in detector /. Although we are dealing 
with Gaussian noise, i.e., h(t; f) is normally distributed, 
the associated PDF for SNR(i; /, Cl) is not expected to be 
normally distributed. It is more peaked than a normal 
distribution and it has broader tails (see Fig. [T|). 

Fig.[T]shows a comparison of SNR(i; /, Cl) for data and 
Monte Carlo. The data corresponds to approximately a 
third of a day of at the beginning of LIGO's S5 science 
run using the Hanford HI and Livingston LI detectors. 
We introduce an unphysical time-shift between the two 
data streams to remove all astrophysical content. Addi- 
tional data processing details are described in Apps. IB II 
andlB"2l as the precise shape of the PDF for SNR(i; /, Cl) 
depends crucially on details of how time series data is 
processed. 

The Monte Carlo histogram is scaled by a normaliza- 
tion factor (derived analytically in App. IB 2]) . which takes 
into account data processing not included in our Monte 
Carlo simulation, e.g., coarse-graining. After applying 
this normalization factor, we find that the standard de- 
viation of the data and Monte Carlo distributions agree 
to better than four significant digits. We conclude that 
data and Monte Carlo are in qualitative agreement. Thus 
we expect that the data are well-behaved enough that we 
can use a Gaussian noise model to assign a detection can- 
didate threshold for SNRr, at least for this choice of pixel 
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FIG. 1: Histogram of SNR(t; /, ti) using 52 s x 0.25 Hz pix- 
els comparing S5 data with an unphysical time-shift (blue) 
to Monte Carlo data (red) and a normal distribution with 
mean = and a = 1 (black). Error bars are too small to see. 




f(Hz) 

FIG. 2: GW strain amplitude spectrum due to PNS con- 
vection in an axisymmetric PNS model at a typical galactic 
distance of lOkpc 0. This plot was generated using the data 
simulated in Q| available at [98l |. 



V. PATTERN RECOGNITION 

In this section we showcase the cross-power statistic 
developed in Sec. IIIII using two different implementations 
(the box search and the Radon search) designed to ad- 
dress two different type of astrophysical scenarios (broad- 
band signals and narrowband signals). 



A. Broadband box search 

We demonstrate how a box-shaped set of pixels can be 
used to search for a broadband GW transient source. For 
illustrative purposes, we consider a toy model based on 
protoneutron star (PNS) convection with a spectrum pro- 
duced in an axisymmetric PNS model assuming a non- 
rotating, 15 Mq progenitor Q (see Fig. [3]). We simu- 
late a d = 4.5 kpc source in the direction of ra = 17hrs, 
deel = 30° at 00:00 GMST on top of simulated detec- 
tor noise comparable to the design sensitivity for initial 
LIGO. We calculate the cross-power statistic 3^r utilizing 
a 200 Hz x 16 s box constructed with the H1L1 detector 
network. We use 4 s x 0.25 Hz pixels, and for each pixel we 
use 20 adjacent segments to calculate a(t;f,Cl), (10 on 
each side). We tile the /i-map and record the 3^r within 
each box. We find that the signal can be recovered with 
SNR(i; /, SI) = 8. The results are summarized in Fig. [3J 



B. Radon algorithm 

Radon transforms are regularly used in imaging prob- 
lems in order to identify line-like features in 2D maps [t| . 
This makes the Radon algorithm useful for looking for 
narrowband GW tracks in ft maps, (see, e.g., [99[). By 



converting from the coordinates (t; /) to impact param- 
eter b and angle 9, a line-like cluster in /i-space is con- 
verted to a peak in Radon space. Thus, this algorithm 
provides a convenient way to search for GW signals that 
manifest themselves as line-like tracks in /t-space. 

For continuous variables, the Radon transform of some 
function g(t;f) is defined as 

TZ(b,9)\g(t;f)} = JdfJ dtg(t;f)6(b-tcos9-f sm9). 

(5.1) 

For discrete variables, the Radon transform becomes 

rc(M)[ff(t;/)]=5>$ff (<;/)■ (5-2) 
t;f 

The weight factors wf. j describe how close a line, parame- 
terized by (b, 9), passes to the center of each /t-pixe l. W e 
use a (modified) Radon transform algorithm from jl0fj |. 
which is one of many possible implementations of the 
discrete Radon transform. 

The estimator for the cross power in a pixel set de- 
scribed by (&, 9) can be written entirely in terms of Radon 
transforms: 



y(n,b,t 



n[Y(t;f,n)&(t;f,ny 
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(5.3) 



The associated variance is 

£ t;/ «/) 2 ^;/,sir 



(n[a(t-j,ny 



(5.4) 



We now consider a toy model of torus excitations from 
long GRBs which are expected to produce line- like 
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FIG. 3: Injection recovery with the box-search algorithm. Top-left: an /t-map of Y(t; f, f2). The injected signal (not visible by 
eye) is indicated with a green arrow. Top-right: an /t-map of SNR(i; /, fi). The injected signal (not visible by eye) is indicated 
with a black arrow. Bottom-left: a histogram of SNRr for a 200 Hz x 12 s box. The blue dashed corresponds to the injection. 
Though the signal is weak in each pixel, the signal obtained by combining every pixel in V is large. Bottom-right: j^r as a 
function of time. 



clusters in /t-space with durations of 2 ~ 200 s. Since we 
are dealing with an elliptically polarized source, y (i; /, fi) 
also depends on inclination angle t and polarization angle 
tp, (see App. IA 5|) . For the sake of simplicity, however, 
we use an unpolarizcd filter, which has been shown to 
do a rea sona bly good job recovering elliptically polarized 
sources jl0l| . We simulate an elliptically polarized wave- 
form (see Tab. on top of simulated detector noise 
comparable to design sensitivity for initial LIGO. Once 
again, we use 4s x 0.25 Hz pixels, and for each pixel we 
use 18 adjacent segments to calculate a(t; /, Cl) (9 from 
each side). The resulting maps of Y(t;f,Cl) and SNRr 
are shown in Fig. 01 For this toy-model, the signal was 
recovered with SNR r = 8.1 at d= 1.7 Mpc. 



parameter 


value 


strain amplitude 


h me = 1.2 x 10" 21 


search duration 


100 s 


(ra, dec) 


(17hr,30°) 


M) 


(0°,0°) 


injection start time 00:00 GMST 


injection duration 


40 s 


distance to GRB 


1.7 Mpc 


(fo,df/dt) 


(998 Hz, 0.03 Hz/s) 



TABLE III: Search and injection parameters used for Fig. [4] 



C. Other algorithms 

In this section we have, for illustrative purposes, pre- 
sented two of the many pattern recognition algorithms 
that may be applied to the problem of looking for fea- 
tures in /t-maps of cross-power. There is a diverse and 
extensive literature devoted to the study of cluster iden- 



tification, (see, e.g., [102L Il03l |). In the next section we 
apply the Radon algorithm to non-GW channels in or- 
der to look for environmental noise artifacts that are 
qualitatively similar to our long-GRB toy model. For 
comparison, we also make use of locust and Hough al- 
gorithms [lfj , which have been proposed as a method of 
identifying long-GRB events in GW data. 
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FIG. 4: Top: an /i-map of SNR(t; /, fi). The long- GRB 
track is a barely discernible diagonal line. Bottom: a map 
of SNRr(f2, b, 0). The signal shows up as a hot spot at 
(6,b) = (39°, -8) with SNR r = 8.1. The sinusoidal patterns 
are (expected) covariances between Radon-map pixels. 



VI. APPLICATION TO ENVIRONMENTAL 
NOISE IDENTIFICATION 

A. Environmental noise in GW interferometers 

While our discussion until now has been focused on the 
detection of GW transients, the same formalism can be 
applied to look for structure in /i-maps of cross-power 
between any two data channels. In particular, it is il- 
luminating to study the cross-power between an inter- 
ferometer's GW-strain channel, (which we denote sqw) 
and a physical environmental monitoring (PEM) chan- 
nel such as a seismometer or a magnetometer channel 
located near the interferometer. Since PEM channels are 
not sensitive to GWs, statistically significant features in 
an ft-map of PEM-sqw cross-power are likely due to 
environmentally-induced noise artifacts. 



Transient artifacts are called "glitches" whereas per- 
sistent narrowband features are often called "lines" or 
"wandering lines" when the frequency slowly changes 
over time. Glitches and wandering lines can be prob- 
lematic for searches for bursts / compact bin ary coales - 
cences and for pulsars respectively, see, e.g., [104 - |107 [. 
(They also produce non-Gaussian noise for our cross- 
power statistic.) It is thus desirable to identify and when 
possible mitigate these noise features. 

In this section we show how the formalism we have 
developed to search for long GW transients can also be 
used to identify glitches and wandering lines in PEM- 
sgw cross-power maps. There are two points we hope to 
make with this digression. First, we shall see that PEM- 
s gw /i-maps are useful for identifying, characterizing 
(and in some cases eliminating) environmentally-induced 
noise. Second, we show that some environmentally- 
induced noise in PEM-sqw /i-maps is qualitatively sim- 
ilar to the GW-transient signature in cross-power maps 
between two sqw channels. Thus, these noise artifacts 
provide a convenient dataset to demonstrate our search 
algorithms with environmental noise events. 



B. Environmental channels at LIGO and Virgo 

In order to facilitate the detection of transient GWs, 
it is necessary to monitor and characterize glitches and 
lines. Efforts to identify and document noise artifacts 
are a major task of the Detector Characterization and 
Glitc h groups wit hin the LIGO Scientific Collabora- 
tion [Tol, [Tol [Toil and the Virgo Collaboration. To as- 
sist in this effort, each LIGO/Virgo observatory is supple- 
mented with hundreds of sensors that monitor the local 
environment. (For an oy ervie w of the LIGO and Virgo 
interferometers, see jll0| and jlllj respectively.) 

Accelerometers measure vibrations such as the motion 
of the beam tubes and of the optical tables that house 
photodiodes; microphones monitor acoustic noise at crit- 
ical locations; magnetometers monitor magnetic fields 
that could couple to the detector; radio receivers monitor 
radio frequency power around the laser modulation fre- 
quencies; and voltage line monitors record fluctuations in 
the AC power. The PEMs are placed at strategic loca- 
tions around the observatory, especially near the corner 
and ends of the interferometer where important laser, op- 
tical and suspension systems reside in addition to the test 
masses themselves. 

For illustrative purposes, we consider a special class 
of noise artifacts induced by passing airplanes. These 
"airplane events," have attractive properties for our pur- 
poses. First, airplane events are relatively well under- 
stood. The existing LIGO airplane veto system (called 
planemon) has been shown to flag airplanes observed in 
microphone channels, and these flags have been shown 
to agree with airplane flight data [l!2l Ill3l | , though the 
existing planemon algorithm does not determine if the 
passing airplane affects sgw(^)- Since we already under- 
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stand a lot about airplane events, it is straightforward 
to assess if our algorithms are consistent with what we 
already know. 

Second, we shall see that airplanes produce a slightly 
curved narrowband /i-map track lasting tens of seconds. 
These tracks are qualitatively similar to GW transient 
scenarios such as van Putten's long-GRBs model, and 
thus they provide an opportunity for us to demonstrate 
our search algorithms on a distribution of signals that re- 
semble unmodcled GW transients; we have a qualitative 
picture of the signal, but it is impractical to model the 
space of all possible signals with a matched filter tem- 
plate bank. 



C. Airplane noise identification 

In order to identify airplane events, /i-maps correlat- 
ing the GW channel with acoustic channels are computed 
in 400 s blocks. Then we take the absolute value of the 
SNR(i; /) map. This step is not necessary for GW studies 
because GW transients can only produce positive-definite 
SNR(t; /) once the phase delay between two sites has 
been taken into account. Transient noise artifacts, on the 
other hand, can produce complex (not positive-definite) 
SNR(t; /) since the cross-power phase depends on the 
coupling of the environmental noise into sgw- By tak- 
ing the absolute value of SNR(i; /), the PDF of SNR(f; /) 
changes from the description in Subsec llVl so we estimate 
it semi-empirically with the assistance of simulation (see 
below) . 

We create a Radon transform of each /i-map of 
|SNR(t;/)|, an example of which can be seen in Fig. [5] 
We are presently interested in qualitative aspects of the 
airplane track, and so we make the simplifying assump- 
tion that a(t; /) = const. Although airplane tracks are 
slightly curved, the approximation of the tracks as lines is 
suitable for a simple identification. We record the bright- 
est spot, i.e., the maximum SNRr, on each Radon map. 
The maps are then ordered according to their maximum 
SNRr- The /t-maps with the highest SNRr are checked 
by eye for airplane tracks. In Fig. [6j we show a histogram 
of SNRr. The red entries are for /i-maps with unam- 
biguous airplane tracks. The green entries are possible 
airplane tracks, yet a visual inspection of these events 
was inconclusive. All other entries (due to non-airplane 
background) are blue. Using Fig. El we set a threshold of 
SNRr > 16 for an event to be tagged as an unambiguous 
airplane. 

In Tab. IIV1 we compare the Radon triggers to the 
planemon triggers. During a span of 3 days, the Radon 
algorithm found 37 unambiguous airplanes, 26 of which 
were also seen by planemon. The 11 events that were 
not confirmed by planemon were confirmed by eye as 
airplane-like. It is possible that they failed to create a 
planemon trigger because the airplane flight path did not 
trigger enough microphone channels. During this same 
period, planemon found 54 unique events, 28 of which 





planemon- flag 


ged not planemon-flagged 


Radon-flagged 


26 


11 


not Radon-flag; 


ged 28 


583 



TABLE IV: Number of airplane event triggers identified by 
the Radon and planemon algorithms out of 648 total ft-mwps. 



were not found by the Radon algorithm as the tracks 
were not coherent in sgw- 

We conclude that the Radon flags appear to be roughly 
consistent with the planemon flags. However, the Radon 
flag differs in two useful ways. First, Radon flags require 
coherence between microphone channels and sqw> an d 
so our algorithm only flags airplanes that contaminate 
SGw(t)- Second, the Radon flag does not need confirma- 
tion between more than one microphone channel so long 
as the signal from one is coherent with sgw- 

For each map deemed to contain an airplane, the in- 
verse Radon transform is performed on the brightest 
spot, and the airplane track is reconstructed as an ft- 
map line as in Fig. [5j The start and stop times of the 
airplane noise in this frequency band are estimated to be 
the times at which the reconstructed track intersects the 
edge of the /i-map. This routine was run on the LIGO 
and Virgo data during recent science runs, correlating the 
GW channel with about 10 microphones at each detec- 
tor. The algorithm identified ~ 10 — 15 airplane events 
at an observatory each day. 

PEM-sgw /i-maps created with real data contain ad- 
ditional environmental noise artifacts besides airplanes. 
In order to estimate the FAR for airplane tracks in "ide- 
alized" noise (where no other environmental artifacts are 
present), we perform Monte Carlo pseudo-experiments 
by scrambling /t-map pixels so as to wash out clustering 
while preserving the empirically observed distributions 
of Y(t;f) and a(t;f). Using the weighted Radon algo- 
rithm, (where a(t; f) ^ const), we estimate that, on av- 
erage, 0.4 events/month are falsely identified as airplanes 
in idealized noise. 



D. Other algorithms 

Having demonstrated how the Radon algorithm can 
recover airplane events in PEM-sgw /i-maps, we now 
demonstrate two additional pattern recognition algo- 
rithms. Our point is to convey the wide array of tools 
available to solve the problem of pattern recognition in 
/t-maps of cross-power. There are typically both advan- 
tages and disadvantages associated with each algorithm, 
which means that each one lends itself to different ap- 
plications, though extensive discussion of the merits of 
different techniques is beyond our current scope. 

In particular, we consider the locust and Hough al- 
gorithms from [Tol | . both of which arc well-suited for 
narrow-band sources. The Hough algorithm is similar to 
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FIG. 5: Top-left: a 400 s-duration map of SNR(t; /) (see Eq. I3.11[l created with 4s x 0.25 Hz pixels and using sgw cross- 
correlated with a microphone. The slightly curved track on the right side of the plot is caused by the Doppler-shifted acoustic 
signal from a passing airplane. Top-right: the associated Radon map. Note the bright spot on the mid-right corresponding 
to the airplane track. Bottom-left: ft-mwp of the reconstructed track using the maximum SNR(t; /) pixel in Radon space. 
Bottom-right: /t-map of the magnitude of SNR(t; /) including a black line corresponding to the veto window. These data are 
from the beginning of LIGO's S5 science run. 



the Radon algorithm, except that it can be extended to 
fit tracks described by arbitrarily high-order polynomials. 
By introducing additional fit parameters, the tracks tend 
to be reconstructed more accurately. However, by adding 
more parameters, the significance of a line-like event with 
little or no curvature can be less than the value obtained 
by the Radon algorithm. 

The locust algorithm is a local wandering algorithm, 
which integrates the /i-map along a chain of local max- 
ima. This algorithm has the advantage that it can re- 
construct arbitrary-shaped tracks without large numbers 
of free parameters. Since it relies on local maxima, how- 
ever, the Radon and Hough algorithms are more robust 
if the GW power is spread diffusely over many pixels. 
Both the locust and Hough algorithms produce a statis- 



tic, which is the integral of cross-power along a track. We 
estimate significance by performing Monte Carlo pseudo- 
experiments in which we randomly scramble the /t-map 
pixels. 



Applying the locust and Hough algorithms to an un- 
ambiguous airplane event, we obtain the reconstruction 
plots shown in Fig. [7J We determine that both the locust 
and Hough algorithms detect the event with a FAR no 
more than 0.04% per 400 s map in idealized noise. 
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FIG. 6: A histogram corresponding to SNRr for above- 
threshold airplane events (in red, dotted outline), near- 
threshold events (green dashed outline) and below-threshold 
background events (in blue, solid outline) in three days of 
LIGO Hanford Observatory data. The dotted black line cor- 
responds to the chosen threshold of 16. The red, green and 
blue distributions are separately normalized to unity for the 
purpose of plotting. Entries with SNRr = record the maps 
which were excluded from this study for failing a "glitch cut." 
Entries with SNRr = 40 record the maps with SNRr > 40. 



VII. COMPARISON WITH OTHER 
TECHNIQUES 

In this section, we compare the proposed excess cross- 
power statistic to matched filtering and to the general 
excess-power statistic from @. It is impractical to ac- 
tually carry out a matched filtering analysis since we do 
not have a template bank from which we can construct 
an arbitrary transient signal. Nevertheless, it is possible 
to perform analytical estimates. 

We consider a signal space characterized by an ft- 
volumc V spanned by N e g independent matched filtering 
templates. The /i-volume is simply the number of pixels 
in our pixel set T. We endeavor to address the following 
question: given a false-alarm probability (FAP) and a 
false-dismissal probability (FDP), what is the minimum 
signal amplitude detectable by either method? Follow- 
ing P, we respectively define thresholds A^ n , A^ n and 
A^f n as the minimum detectable amplitudes by our cross- 
power (CP) search, by the optimal total excess-power 
(EP) search Q, (which includes both cross-power and 
auto-power terms) and by a matched filter search (MF). 
The ratio of the corresponding amplitudes is (by defini- 
tion) the efficiency of the excess power statistic compared 
to matched filtering: 



?7epmf(FAP, FDP, N e g, V) = A™?/A^ n (7.1) 
,7cPEp(FAP,FDP,y) = AZjAgl. (7.2) 

To calculate these thresholds, signal and noise distri- 
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FIG. 7: Upper panel: an airplane track in a map of 
SNR(t; /, Q). Middle panel: the reconstructed track as de- 
termined by the locust algorithm. Bottom panel: the recon- 
structed track as determined by the Hough algorithm using a 
second-order polynomial. 



butions for CP and EP are generated using Monte Carlo 
simulations. Throughout this section we assume station- 
ary Gaussian white noise; simulated signals are char- 
acterized only by their amplitude, polarization and ft- 
volumc. Other characteristics such as frequency content, 
evolution with time, etc. are not relevant for this white- 
noise calculation. Following [1] , we approximate the MF 
threshold as 



.4 



MF 



A^ n (FAP/JV eff) FDP,l/2). 



(7.3) 



In Fig. [51 we plot ?7epmf as a function of V and N e g. We 
see that ??epmf ^ 50% over the range of parameter space 
considered. 

To compare the CP method to the EP method, 
we calculate ?7cpep for the Hanford-Livingston and 
Hanford- Virgo networks averaging over an isotropically 
distributed population of unpolarized GW sources. In 
Fig. [9] we plot 7?cpep as a function of ft- volume V. The 
CP technique is highly efficient (??cpep ~ 77%) even at 
small values of V, but it becomes increasingly more effi- 
cient at higher V. 
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FIG. 8: The relative detection efficiency ?7epmf of excess- 
power compared to matched filtering as a function of ft- 
volume V and the number of effective templates N e a (for 
FAP = 1CT 2 and FDP = 1CT 1 ). We assume Gaussian white 
noise. For a given value of V, there is maximum number 
for N e g above which the templates will not be independent. 
Therefore, at high values of N e g and small values of V, there 
is an unphysical region where rj > 1. A realistic excess-power 
(and hence cross-power) search might have V ~ 0(100), 
which corresponds to an efficiency of ~ 50% compared to 
matched filtering. 



Ir 
0.95 
0.9 
0.85 



0.8 
0.75 

0.7 

0.65 

0.6 . 
10 




10 10 
Volume (V) 



a novel framework, which can be used to look for GW 
transients on timescales of seconds to weeks. This frame- 
work, which is a generalization of the GW radiometer [tJ , 
utilizes /i-maps of GW-strain cross-power using two or 
more spatially separated interferometers in order to look 
for statistically significant clustering. A comparison of 
simulated detector noise with time-shifted data revealed 
that /i-cross-power-maps made with real interferometer 
data are well-behaved (for at least one pixel size) sug- 
gesting that the threshold for candidate events can be 
determined analytically. 

We illustrated how different pattern-recognition tech- 
niques can be used to identify GW signatures in ft- 
maps. We demonstrated some of these techniques us- 
ing /i-maps generated using a GW-strain channel cross- 
correlated with a LIGO microphone channel and we pre- 
sented a novel technique for the identification of envi- 
ronmental noise transients in GW interferometers. Wc 
assessed the effectiveness of our proposed statistic com- 
pared to matched filtering and other excess-power strate- 
gies. We found that our strategy is highly effective when 
looking for long GW transients whose precise waveforms 
are not known. 

Many of the sources we considered here arc plausi- 
ble targets for the advanced detector era. If long GW 
transients are, in fact, detected, we shall gain invalu- 
able information about objects and processes, for which 
we currently possess only preliminary models, e.g., long 
gamma-ray bursts. If, on the other hand, no long GW 
transients are detected, we expect that some models pre- 
dicting relatively large strain amplitudes (e.g., ||) may 
be ruled out or constrained. Third-gener atio n detectors 
such as the proposed Einstein Telescope |l!4| can apply 
our long GW transient algorithm to probe still fainter 
sources. 



FIG. 9: The relative detection efficiency comparing the gen- 
eral excess power and cross-power methods for FAP = 10~ 2 
and FDP = 0.1 in simulated Gaussian white noise. The red 
circles are for the HI VI detector pair and the blue squares 
are for the more nearly aligned H1L1 detector pair. The er- 
ror bars represent one standard deviation using 100 trials. 
The cross-power method has a relative detection efficiency of 
at least 70% with respect to the excess power statistic and 
the efficiency increases with V. 



VIII. CONCLUSIONS 

We reviewed a variety of compelling scenarios for long- 
GW transients including protoncutron star convection, 
accretion disk fragmentation/excitations, rotational in- 
stabilities in neutron stars, r-modes, pulsar glitches and 
soft gamma repeater flares. Many of the models we 
considered predict strain amplitudes detectable in the 
advanced-detector era. 

Next, inspired by stochastic analyses, we introduced 
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Appendix A: Derivations 

We describe how a GW source can be characterized by 
its power spectrum H(t;f), we construct an estimator 
Y(t; /) for H (t; f) and calculate the associated variance. 
We construct an estimator for the variance. Finally, we 
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construct the filter function for an clliptically polarized 
source. 



1. Introduction and notation 

Working in the transverse-traceless gauge we write 
down the general form of a GW field, which can depend 
on direction f2, polarization state A and frequency / (see, 
e.g., HSU!): 

h ab (t,x) = V f°° df [ dCiet b (Ci)h A (f,Ci) 

a J-°° J s 2 (Al) 

2irif(t+Q-x/c) 

Here x and Cl are defined in the reference frame with the 
origin fixed at the center of Earth, but not rotating with 
the Earth. The indices ab run over a Cartesian coordinate 
system. We define unit vectors 

9 = cos 9 cos 4> x + cos 9 sin cf> y — sin 9 z (A2) 
4> = — sin x + cos <f> y (A3) 
Cl = sin 9 cos <f> x + sin 9 sin <fiy + cos 9 z ( A4) 

such that {9, <j>, Cl} form a right-handed coordinate sys- 
tem and the rotational axis of the Earth points along z. 
The two GW polarization tensors can be written as (see, 
e.g., 0, HI 123): 



e+ (fi) = 9 <g> 6-<j>® cj) 
e ab(^) = ® $ + 4) ® e. 



(A5) 
(A6) 



Since we are looking for GW transients, we restrict 
our attention to point sources for which /ia(/, Cl) = 
h A (f)S(Cl — Cl ). We perform the integral over Cl and 
obtain: 



h ab (t,x)=J2 dfe£ b (Cl )h A (f)e 

A J -°° 



2nif(t+Cl -x/c) 

(A7) 



For simplicity, we henceforth replace fio with Q. It fol- 
lows that the GW strain in detector I is given by 



hi(t) = Y, dfh A (f,n)e 

A J - aa 



eMdf(t) 
(A8) 

where dj b (t) is the detector response tensor at time t: 

di(t) = i (Xi(t) g) Xj(t) - y/(t) ® F/(0) ■ (A9) 

Here, the two detector arms lie along the X(t) and F(i) 
axes, which are time-dependent due to the rotation of the 
Earth. 



We now consider a finite stretch of h] (t) and take the 
discrete Fourier transform of Eq lA8l to obtain: 



h I {t;f)=Y,hA(t\f,Cl)e 



2itiftt-x I /c rpA 



Fj A {t;n). (A10) 



where we define the "antenna factors" (see e.g., (921) to 
be 

Ff(t;ti) = e&(ti)df(t). (All) 
We define the GW strain power spectrum to be 

(h A (t;f)h A ,(t;f)) = ±H AA ,(fJ), (A12) 

where the factor 1/2 comes from the fact that H A A'{t'i f) 
is the one-sided power spectrum. Here we use the discrete 
Fourier transform defined in Eg. 13.21 (see also Tab. HI)) . 

It is convenient to characterize the source with a single 
spectrum that includes contributions from both + and x 
polarizations. We therefore define 



H(t;f) = Tr [HAA,(t;f)] 



(A13) 



so as to be invariant under change of polarization bases. 
This definition is a generalization of the one-sided power 
spectrum for unpolarized sources found in 0,112, Our 
goal now is to derive an estimator for H(t] f) in a data 
segment over which it is presumed to be constant. 



2. Derivation of Y 

Let s/(t) = hj(t) + ni(t) be the strain time series from 
detector /, where hj{t) is the GW strain and nj(t) is 
the detector noise. Following 0, H2, HH, we combine 
the strain time series from two spatially separated detec- 
tors, Si(t), Sj(t), to construct an estimator for GW-powcr 
H(t] f) for a point source at a sky position 17, 

Y(t;f,Cl) = 2 Re f, A)SJ(t; /)Sj(t; /)] (A14) 

where Qu(t;f,Cl) is some filter function to be deter- 
mined below. We take the real part to ensure physicality 
of the estimator. The expectation value of Y(t; /, Ci) is 
given by 

(Y(t; /, f2)) = 2 Re [Q„(t; /, /) hj(t; /)) 

(A15) 

since, by assumption, there is no correlation between sig- 
nal and noise and also no correlation between noise in 
two spatially separated detectors. 

Combining Eqs. IA101 IA15I and IA12I we get 

(Y(t; /, f2)) =2 Re \Qu(t; /, Cl) £ \h aa , (t; /) 

aa> 1 



e -2^/(a.(S 7 -x J )/c) pA^. F A>( f . 



(A16) 
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In order to simplify the form of HAA'(t; f) we now con- 
sider unpolarized sources, for which 



H AA ,(t;f) = ^H(t;f)5 AA , 



(A17) 



(We consider the case of polarized sources in Subscc. lA 51 ) 
For unpolarized sources, 



(Y(t;f,n)) = lRe Q„(t; H(t; f) , 



5^(t;n)F*(t;ft) 



where we have defined Axu = xi — xj. 

We desire that (Y) = i? (i; /), which implies: 



-2mfn-Axij/c 



(IJ 



(f;/,n) 



2 g27Tl/0- Axjj/c 

£ A F/(t;n)F/(t;fi)' 



(A18) 



(A19) 



By setting Quit] /, 17) thusly, we account for the phase 
difference between detectors / and J ensuring that the 
bracketed quantity in Eq. IA18I is real. We also account 
for the detector pair efficiency. 

Finally, we define (unpolarized) pair efficiency as 



Z A 

which enables us to rewrite the filter function as 

1 



Q LI (t-j,n) = 



= 2iri/n-Ax/j/c 



(A20) 



(A21) 



Since Y(t; f, SI) cx Q(t; /, 0) and Q cx l/e/jfi; /, CI), it 
follows that Y(t; /, f2) cx l/e/j(t; SI). This can be under- 
stood as follows. If we observe a modest value of strain 
power from a direction associated with low efficiency, we 
may infer (if the signal is statistically significant) that 
the true source power is much higher because the net- 
work only "sees" some fraction of the true GW power. 



3. Variance of the estimator 

We derive an expression for the variance of Y{t; /, f2), 
avitj^) 2 ee (Y(t-J,{l) 2 ) - (Y(t;f,Q)) 2 . In searches 
for persistent stochastic GWs, the second term is usually 
omitted and the first term is simplified by assuming that 
signal in each pixel is small compared to the noise. Such 
small signals arc extracted by averaging over a very large 
number of segments (see, e.g., [Hi). Since we are dealing 
with transients, however, the signal may be comparable 
to the noise and so we can not neglect any terms in our 
calculation of a\. 

To begin we define a new (complex-valued) estimator 
that will be handy in our derivation of a\: 



Our GW power estimator Y(t; /, SI) is simply the real 
part of W{t;f,n): 

Y(t; f, fi) = i (W(t; /, Si) + W(t; /, SI)*) . (A23) 

For notational compactness, we shall omit the arguments 
of W(t; f, SI) in the remainder of this derivation. It fol- 
lows that the variance of Y(t; f, SI) can be written as 



n 2 1 



(w 2 ) - (wy + (to - (w*y + 



2a w 



where 



oiy = (\W\')-\(W)\' 



(A24) 



(A25) 



Now we evaluate the three terms in Eq. I A24I beginning 
with aly. We obtain 



{si(t;f)sj(t;f)Si(t;f)3j(t;f)) 



(5/(*;/)5j(*;/))<Sj(t; /)«}(*;/)) 

2 



(*;/,«) 



(A26) 



For mean-zero Gaussian random variables, we can ex- 
pand the four-point correlation into a sum of products of 
two-point correlations. We substitute s = h + n and set 
signal-noise cross terms to zero along with noise-noise 
cross terms from different detectors. The variance be- 
comes 

a 2 w (t; f, SI) =4 [(h*j(t; /)^(t; /)) (hj(t; f)h*j(t; /)) + 

(&K*5 /)&/(*;/)) (nj(t;f)n*j(t;f)) + 
(hj(t;f)h*j(t;f)) (n* I (t;f)n I (t;f))+ 
(n}(t;f)nj(t-J)) (nj(t; /)«*(<; /)) 



(A27) 



Evaluating the four terms in Eq. IA271 we obtain 



<7^(t;/,fi) 



e n {t;n)ejj(t;n)H(fjy 



H 



(t; /) (eu(t; Cl)Nj(t; f) + ejj(t; ft)JV}(f; /)) (A28) 



+ N I (t;f,)N J (t;f) 



Quit; f, to) 



where e is defined in Eq. IA20I and where Nj(t; f) is the 
one-sided noise-power spectra: 



W(t; f, n) ee 2 Q u (t; f, n)Sj(t; f)ij(t; /). (A22) 



Nj(t;f) = 2\ni(t;f)\' 



(A29) 
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Using the same line of reasoning, we calculate the re- 
maining terms in Eq. IA24I 

(W 2 ) - (W) 2 = (W* 2 ) - (W*) 2 = H(t; f) 2 . (A30) 

Combining Eqs. IA241 and lA"30l we conclude that 



a Y = -[a 2 w + H(t;f) 2 }. 



(A31) 



The factor of 1/2 comes about from the fact that 
Y(t; f, fi) is real whereas W(t; f, ft) is complex. We note 
that in the small-signal limit H(f) — > and the variance 
reduces to the canonical stochastic result 1921]: 



N T (t;f)Nj{t;f) Qu(t;f,ty , (A32) 



1 



4. Expectation value of a Y 

Our estimator for the variance of Y is given by 

1 _ .2 



#(t;/,n) = - Qu(t;f,n) Pf\f)PT(f), (A33) 



,ad j(/)p adj ( 

where Pj is the average auto-power in neighboring pixels: 



p; dj (/)^2|5 / (/)r 



(A34) 



The ovcrline denotes an average over neighboring pixels. 
By averaging over neighboring pixels, we assume that the 
detector noise in any given pixel can be characterized by 
looking at its neighbors. This assumption is discussed 
below. 

Now we calculate the expectation value of our estima- 
tor for variance dy given in Eq. IA33l in order to compan 



it to the theoretical variance given in Eqs. IA311 and lA"28l 
Eqs. 13.91 and 13.101 together imply 



(5f(*;/,A)> =2 Qu(t;f,n) 



(A35) 



Using Equation [A26] to write the expectation value of a Y 
in terms of the theoretical value of ct 2 ^ , we find 



<#(*;/, ft)) =- 



<7^(i;/,n) 



(*;/,«) 



(sf *(t; f)Sf (i; f))(Sf(t; f)sf*(t; /)) 



a 2 w (t;f,n)+ (W) 



a 2 w {t;f, Cl) + H(t;ff 



(A36) 



Since this is the theoretical variance from I A3 11 we con- 
clude that (fry) = a Y . Thus, Eq. 13.91 provides an unbi- 
ased estimator for a 2 ,. Here we have assumed that the 



noise and signal are comparable in neighboring segments. 
This assumption can fail for rapidly changing, high-SNR 
signals and also for highly non-stationary noise, and so 
additional work may be required to estimate a in these 
situations. 



5. Elliptically polarized sources 

A variety of long-transient GW sources are expected 
to be ellipticall y po larized (e.g., long GRBs d, Q and 
pulsar glitches jll5l |). Elliptically polarized sources are 
parameterized by two angles. The inclination angle i is 
the angle between the rotational axis of the source and 
the observer's line of sight and the polarization angle 
ip describes the orientation of the rotational axis in th e 
plane perpen dicu lar to the line of sight (see, e.g., |116| ). 

Following |117| . we characterize an elliptically polar- 
ized source with the so-called canonical amplitudes: 



A, 



where 



A x 



( A + cos2V> \ 
A + sin 2-0 
-A x sin 2^ 
\ A x cos2V>, / 



{ha/2) (1 + cos 2 l) 

ho COS L 



(A37) 



(A38) 
(A39) 



We have set the initial phase </>o = for the sake of sim- 
plicity. (Ultimately, we are concerned with the average 
cross-power over many cycles and so the initial phase is 
unimportant.) Here ho is the stra in amplitude. 
Next, we define the tensor [l!7l | 



K b (t) 



e* 6 cos [tj>(t)] 
e+ sin [<£(*)] 
Ve a >n[<Hi)] J 



(A40) 



whe re 6 {t) describes the phase evolution of th e sig- 
nal [l2ll | . Now we write the GW equation as |ll7l ] 



hab(t) 



(A+ cos(2^) cos [<j>(t)] - A x am(2ip) sin [4>(t)}) e+ + 

(A+ sin(2V) cos [<f>(t)] + A x cos(20) sin [(f>(t)]) e* b . 

(A41) 

We Fourier transform a finite stretch of GW signal to 
obtain the coefficients h,A(f) (see Eq. IA1|) : 

~ h+ (t-f) = [A+ cos(2V;) + iA x sin(2VQ] 6 ffo 



hx(t;f) 



2N S 

[A+ sin(20) - iA x cos(20)] 5f fo 
2N S 



(A43) 
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We have ignored negative frequencies, which will play 
no part in our subsequent calculation of one-sided GW 
power. Also, we have expanded <p(t) in a Taylor series, 



ht + \ft 2 ) +0(t 3 ). 



(A44) 



We assume that / is approximately zero over the segment 



J 



duration, T. In principle, this formulation could be ex- 
tended to broadband sources, but we expect elliptically 
polarized sources (associated with spinning objects and 
a binary objects) to be narrowband. 



Combining Eqs. IA12I and IA42| we calculate the one- 
sided power spectrum: 



H AA >(t;f) 



°ffo 

"2JV? 



A\ cos(2i/;) 2 + A\ sin(2V>) 2 {A\ - A\ ) cos(2V>) sin(2^) 
{A\ - A\ ) cos(2V>) sin(2V) A\ sin(2^) 2 + A\ cos(2V>) 2 



A+A, 



-A + A> 




(A45) 



The imaginary off-diagonal terms corresponds to the structed from the cross power spectrum of two GW- 



phase delay between + and x polarization states. 
Combining Eqs. I A13l and IA"45| we obtain 



H(t; f) 



2N? 



(Al + A 2 x )S fh . 



(A46) 



Our goal, once again, is to find an estimator for H(t] /). 
We assume that the estimator for H(t; /) can be con- 



strain channels multiplied by an appropriate filter func- 
tion, Qu(t; f,Cl,i,tp): 

Y(t; /, O, 4, VO = 2 Re [gj(t; f)sj(t; f) Q u (t; /, ft, i, V>) . 

(A47) 

Plugging Eq. lA"45l into Eq. IATgI we obtain 



2 Re 



{FtFj [A\ cos(2^) 2 + A\ sin(2^) 2 ] 
F+F* [(A 2 + -A 2 x )cos(2i})sm(2ij)~iA + A x ] + F? Fj [(A\ - A\ ) cos(2t/>) sin(2^) + iA + A x ] 
F?F$ [A\ sin(2V.) 2 + A 2 X cos(2^) 2 ] ) e -W(ft-A*Wc) g Jj(t . /; 0) t) ^) , 



(A48) 



r 



We write the polarized filter function as: 



!IJ 



(t; /, ft, 4, V-) - e^iW-ASrj/c+ir, / e fa ^ ^ ( A49 ) 



where 77 is an angle arising from from the phase de- 
lay between + and x polarizations. By requiring that 
(Y(t; /, ft, l, VO) = H(t; /), it follows that 



e(Ct, l, ip) = F+Fj[a 2 + cos(2^) 2 + a 2 sin(2^) 2 ]+ 

{FfFj + F?F}){a 2 + - a 2 ) cos(2?/>) sin(2</>)+ 
FfFj {a\ sin(2?/>) 2 + a 2 cos(2V>) 2 ) 
+ ia+a x (F?F+ -F+F}) 

/(4 + °x). 

(A50) 



and 

7?(Q,4,V) 



-phase I 



F/FJ" [< cos(2V>) + a 2 sin(2r/0 ]+ 



{FfFj + FfF+){a 2 + - a 2 ) cos(2^) sin(2^)- 
F*F* {a 2 + sin(2iA) 2 + a 2 cos(2^) 2 ) 



where we have defined 



a+ = A + //i 
a x = A x /ho. 



(A51) 

(A52) 
(A53) 



As a sanity check, we note that an elliptically polar- 
ized source with pure +-polarization yields a sensible ef- 
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ficiency 



e(n,t = 90°, tp = 0) 



(A54) 



Appendix B: Data processing details 

1. Overview 

In this we describe the detailed procedure for the cre- 
ation of the 52 s x 0.25 Hz coarse-grained data in Fig. [T] 

• Since n(t) is drawn from a zero-mean, unit- variance 
normal distribution, the frequency domain data 
will be zero-mean with a standard deviation scaled 
by yf N/2 for the real and imaginary parts sepa- 
rately where N is the number of sampling points. 

• To calculate the cross-spectral density, CSD, data 
are first Hann windowed and zero padded in the 
time domain. Discrete Fourier transforms are mul- 
tiplied to form cross-power, then summed (with a 
window that tapers at the ends) to yield a CSD 
with 0.25 Hz resolution. 

• To calculate each Pj(/), data are broken into 4 s 
blocks with 50% overlap. Each block is Hann win- 
dowed, and the resulting power samples from all 
blocks are summed so that Pi(f) also has 0.25 Hz 
resolution. 

• Hann windowing reduces the overall power by a 
factor of 3/8. To compensate both the CSD and 
Pi (/) are corrected for this factor. 

• Without windowing or overlap each power estimate 
from each data block will be x 2 -distributed with 2 
degrees of freedom and scaled by N/2, and so the 
average of M blocks are x 2_ distributed with 2M 
degrees of freedom and scaled by N/2M. 

• Overlap changes the effective number of degrees of 
freedom in a way that can be calculated numerically 
(see App.lB2|). 

• Without windowing, each cross-power product will 
have zero-mean and a standard deviation scaled by 
N/y/2. Summing M values gives each CSD esti- 
mate a standard deviation scaled by Ny/M/2. 

• These scaling behaviors indicate how we could ex- 
pect SNR to behave, except for the non-trivial ef- 
fects of windowing, zero-padding and weighted av- 
eraging. These somewhat subtle effects are evalu- 
ated explicitly in App. 



2. Calculation of the normalization factor 

Here we derive the stretch factor applied to the signal 
to noise ratio (SNR) histogram in Fig. Q] generated using 



simulated Gaussian noise. Let {x n } and {y m } be real, 
discrete, time sequences corresponding to two indepen- 
dent data streams with each x n and y m sampled from 
a (0,1)-Gaussian distribution. The time-domain {x n } 
and {y m } can be transformed into frequency-domain se- 
quences of N complex numbers {X p } and {Y q } using the 
discrete Fourier transform (DFT), defined here as: 



N-l 

X p = At ^2 

cxp 

n=0 
N-l 

Y q = At ^ eX P 

m=0 



—2nipn\ 



N 



and 



'2iriqm\ 

-J Vn 



(Bl) 



N 



The inverse DFT can be used to express the original 
discrete-time signal as 



C = AfJ2 exp ( ^rp) x p 



N-l 



Vm 



cxp 



q=0 



V N 
( 2-Kiqm\ 



and 



(B2) 



Y„ 



To see the effects of zero padding, we define the real, 
2N-periodic discrete-time signal x = {x n ,0 n } where n 
is a sequence of N-many zeros. Then 



where 



2JV-1 

A< exp 

n=o 

N-l 

At exp 

n=0 



My 2 *" 1 ext) 

2 l^p=0 LX P 



2iripn\ _ 



2N 



—2 / nipn\ 

I Xr, 



(B3) 



2N 



' 2nipn \ y 
, 2N ) 



¥£p=o~Vl) p exp(^)X p 



with < n < N - 1 



(B4) 



From Eq. IB4I we conclude that the X p with < p < 
2N — 1 are not linearly independent. However, the 
X 2p (0 < p < N - 1) and X 2p+ i (0 < p < N - 1) 
are separately linearly independent. In the calculations 
below, we use only the even terms to express both the 
even and odd X/~ terms. Therefore, taking into account 
the zero-padding, we can write 



N-l 



x n = A/ ^ exp 



p=0 
N-l 



2ixi(2p)n 
27V 



X 



2;j 



= A /E exp(^U p 

p=0 ^ ' 

By comparing x n in Eqs. IB2IIB5I we see that 
X 2 p = X p 



(B5) 



(B6) 
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In addition to the zero padding, the windowing effect can 
be expressed by defining u = {u n ,0 n } = {x n f(n),O n } 
with the discrete Fourier transform given by 



2 \~^ 1 ( —2nikr^ 
U k = At ex P I — — — I 



n=0 
N-l 



At exp 

N-l 

At exp 



n=0 



2N 

—2irikn 
2N 

—2-Kikn 
2N 



(B7) 



f(n)x n 



where x n is given by Eq. IB5I and f(n) is the Hann 
window defined by 



/(») = ~ 2 



1 — cos 



2nn 



a=-l 

with < n < N - 1 



1 /-l\ |a| [2mna\ 



(B8) 



where Q(t; f, Ct) is the filter function, we observe that 
in order to explain the SNR distribution obtained from 
the processed data we need to understand the statistical 
consequences of the code that calculates the power 
spectral density PSD^ (or PSD y ) and the CSTD xy . 

Since the CSD is the Fourier transform of the 
cross-covariance function, in order to characterize the 
numerator of SNR we need to study terms of the form 
X^Yk- Without windowing and without zero-padding 
the X^Yk term depends only on the frequency k as 
the following result shows (obtained using Eq. IB 1 1 and 
Eq.El. 



n=0 m=0 p=0 9=0 ^ 



exp 



-2~Kim(k — q) 
TV 



JV-l N-l 

E E S PkS qk X* p Y q 

p=0 q=0 



(B12) 



After some algebra Eq.ESJ Eq.LBgJ Eq.[BZ]and Eq.[B8] 

give 



A/At 



iV-lJV-l 1 

EEE^ 

n—0 p—0 a— — 1 



exp 



-2irin(k/2 — p — a) 
N 



(B9) 



Similarly, we can express the second data stream by v = 
{v m ,Q m } = {y m f(m),O m } with < m < N - 1. In 
analogy with IB9I the discrete Fourier transform is given 
by 



A/At 



JV-1JV-1 1 

EE En 

rn=0 q=0 b=-l 



exp 



-2nim(k/2- q- b) 
N 



(BIO) 



The sequences {u n } and {v n } are the real, 2N-periodic, 
discrete-time signals that are used in the code for the cal- 
culation of the contribution of the (CSD) to the (SNR). 
From the definition 



SNR(t; /) =Re 



CSD^i;/) 



y/PSD x (t;f)PSD y (t;f) (Bn) 
phase(Q(i;/,n)) , 



However, when we take into account the effects of 
windowing and zero-padding, the U^Vk term, contrary 
to the X^Yk term, has contributions from a wider part 
of the frequency domain according to the following result 



N-1N-1N-1N-1 1 1 



^* = zEEEEE e*;m~ 

n=0 m=0 p=0 q=0 o=-l 6=-l V 



|a| + |6| 



(2mn(k/2 - p- a)\ f -2mm(k/2 - q - b) 

6XP { N J CXP { N 

(B13) 



To imitate the way the code is calculating the CSD, 
we define 



-13 



Zj= E "W+M 



j+k 



k=-13 



where 



= (5, 1,1, ...,1,5) 



(B14) 



(B15) 



25 — times 



It is easy to see that the expectation value of Zj is 
zero, since the X q and Y q are statistically independent 
with mean zero. To evaluate the standard deviation of 
Zj /V26, where the 26 is a normalization constant due to 

averaging, we then need to evaluate (J Z*Zj/2&) (where 
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there is no sum over j). After some calculations, the final 
expression for (Z*Zj) becomes 

k,k' ,a,a' ,£,/3' 



(\a\ + \b\ + \a'\ + \b'\) 



(B16) 



Using N = 212, 992, corresponding to 52 s of data sam- 
pled at 4, 096 Hz, the numerical result turned out to be 



-(Z;Z,)= 0.7208. 



(B18) 



where 



sin[(7r/2)(fc-fc' + 2a-2a')] 
= sin [(7r/2iV)(fc - k' + 2a - 2a')] 



exp 



-i%(a-a') [ 1- — 



(B17) 



Multiplying this by the normalization factor of 8/3 due 
to Hann window we get 1.922 — the normalization factor 
applied to MC in Fig[TJ 
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segments. 

In order to chose a suitable number of neighboring pixels 
to average over, one must typically take into account 
the stationarity of the detector noise. This discussion, 
however, is outside our current scope. 
Here we also assume that the probability density func- 
tion for each pixel is the same, which is to say that the 
noise and signal are approximately stationary. 
The reader is advised to u se ca ution when comparing 
our formulas with those in [ll7j |. Following [7l [92. |93(] . 
we use e^b to denote the "detecto r pol arization basis" 
defined by the ecliptic plane. In [1171 ]. however, this 
quantity is denoted e^t an d e-ab ls the basis in which 
the plus polarization is maximal. 



